function Fisher = ecmnfish(Data, Covar, InvCovar, MatrixFormat)
%ECMNFISH Evaluate Fisher information matrix for ECMNSTD.
%	NUMPARAMS x NUMPARAMS Fisher information matrix based on current maximum
%	likelihood parameter estimates.
%
%   Fisher = ecmnfish(Data, Covar);
%   Fisher = ecmnfish(Data, Covar, InvCovar);
%	Fisher = ecmnfish(Data, Covar, InvCovar, MatrixType);
%
% Inputs:
%   Data - NUMSAMPLES x NUMSERIES matrix of observed multivariate normal data
%		with NaNs to represent missing values.
%   Covar - NUMSAMPLES x NUMSERIES matrix of estimated covariance of Data.
%
% Optional Inputs:
%   InvCovar - Inverse of covariance matrix, i.e., inv(Covar).
%	MatrixFormat - String that identifies which parameters to be included in the
%		Fisher information matrix. The default method is 'full'. The choices
%		are:
%		'full' - (default) Compute the full Fisher information matrix.
%		'meanonly' - Compute only components of the Fisher information matrix
%			associated with the mean estimates.
%
% Outputs:
%   Fisher - NUMPARAMS x NUMPARAMS Fisher information matrix based on current
%		parameter estimates, where NUMPARAMS = NUMSERIES * (NUMSERIES + 3)/2
%		if MatrixFormat = 'full' and NUMPARAMS = NUMSERIES if MatrixFormat =
%		'meanonly'.
%
% WARNING: If calculating the full Fisher information matrix, this routine
%	is VERY slow.
%
% See also ECMNHESS, ECMNSTD, ECMNMLE.

%	Author(s): R.Taylor, 4-21-2005
%	Copyright 2005 The MathWorks, Inc.
%	$Revision: 1.1.6.2 $   $Date: 2005/06/17 20:23:20 $

% Step 1 - check arguments

if nargin < 2
	error('Finance:ecmnfish:MissingInputArg', ...
		'One or more of the required input arguments Data and Covar is missing.');
else
	if isempty(Data)
		error('Finance:ecmnfish:EmptyInputData', ...
			'The required input argument Data is empty.');
	end
	if isempty(Covar)
		error('Finance:ecmnfish:EmptyInputCovar', ...
			'The required input argument Covar is empty.');
	end

	NumSeries = size(Data,2);

	if ~all(size(Covar) == [NumSeries, NumSeries])
		error('Finance:ecmnfish:IncompatibleCovar', ...
			'The covariance matrix Covar has wrong dimensions.');
	end
end
if nargin < 3 || isempty(InvCovar) || ~all(size(InvCovar) == size(Covar))
	InvCovar = inv(Covar);
end
if nargin < 4
	MatrixFormat = 'FULL';
end

% Step 2 - initialization

MatrixFormat = upper(MatrixFormat);
if ~any(strcmp(MatrixFormat,{'MEANONLY','FULL'}))
	warning('Finance:ecmnfish:UnknownFormatString', ...
		'Unknown MatrixFormat string. Will use default FULL.');
	MatrixFormat = 'FULL';
end

if strcmp(MatrixFormat,'MEANONLY')
	NumParams = NumSeries;
else
	NumParams = NumSeries + (NumSeries * (NumSeries + 1))/2;
end

Fisher = zeros(NumParams,NumParams);

% Step 3 - do partials wrt Mean

for i = 1:NumSeries
    for j = 1:i
        Fisher(i,j) = InvCovar(i,j);
        Fisher(j,i) = Fisher(i,j);
    end
end

if strcmp(MatrixFormat,'FULL')

% Step 4 - do partials wrt Covar

	GradC1 = zeros(NumSeries,NumSeries);
	GradC2 = zeros(NumSeries,NumSeries);

	i = NumSeries;
	for i1 = 1:NumSeries
		for j1 = 1:i1
			i = i + 1;

			GradC1(i1,j1) = 1.0;                    % do dC/dtheta(i)
			GradC1(j1,i1) = 1.0;

			j = NumSeries;
			for i2 = 1:NumSeries
				for j2 = 1:i2
					j = j + 1;

					if (j <= i)
						GradC2(i2,j2) = 1.0;        % do dC/dtheta(j)
						GradC2(j2,i2) = 1.0;

						Temp1 = InvCovar*GradC1;
						Temp2 = InvCovar*GradC2;

						Fisher(i,j) = 0.5*trace(Temp1*Temp2);
						Fisher(j,i) = Fisher(i,j);                        

						GradC2(i2,j2) = 0.0;        % undo dC/dtheta(j)
						GradC2(j2,i2) = 0.0;
					end
				end
			end

			GradC1(i1,j1) = 0.0;                    % undo dC/dtheta(i)
			GradC1(j1,i1) = 0.0;
		end
	end
end
